function r=markov(data, min_max, direction, distance)
    vmax = min_max(2);
    vmin = min_max(1);
    nbins = vmax-vmin + 1;
    [M,N] = size(data);
    switch (direction)
        case 'x'
            spacex=distance;
            spacey=0;
            addx=distance;
            addy=0;
        case 'y'
            spacex=0;
            spacey=distance;
            addx=0;
            addy=distance;
        otherwise
            error('unknown direction');
    end
    r = zeros(nbins, nbins);
    for i = 1:M-spacey
        for j = 1:N-spacex
            r(data(i,j)-vmin+1, data(i+addy,j+addx)-vmin+1) = r(data(i,j)-vmin+1, data(i+addy,j+addx)-vmin+1)+1;
        end
    end
    for i = 1:nbins
        s = sum(r(i,:));
        if (s ~= 0)
            r(i,:) = r(i,:) ./s;
        end
    end
end